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We derive and analyze a hierarchy of approximations to the strongly correlated limit of the Hohenberg-Kohn 
functional. These "density representability approximations" are obtained by first noting that in the strongly 
correlated limit, A'^-representability of the pair density reduces to the requirement that the pair density must 
come from a symmetric TV-point density. One then relaxes this requirement to the existence of a representing 
symmetric fc-point density with k < N. The approximate energy can be computed by simulating a fictitious 
A;-electron system. We investigate the approximations by deriving analytically exact results for a 2-site model 
problem, and by incorporating them into a self-consistent Kohn-Sham calculation for small atoms. We find 
that the low order representability conditions already capture the main part of the correlations. 



I. INTRODUCTION 

Kohn-Sham (KS) density functional theory (DFT) is 
currently the most widely used ab initio electronic struc- 
ture model which is applicable to large and complex sys- 
tems ranging from condensed matter over surfaces to 
nanoclusters and biomolecules. With the advent of lin- 
ear scaling algorithms, the key factor limiting the accu- 
racy of predictions is the choice of underlying exchange- 
correlation functionals^^. These functionals model the 
correlation structure of the system as a universal func- 
tional of its underlying one-body density. While highly 
successful in many instances, these functionals exhibit 
known failures, both specific, like incorrect filling order 
and lack of binding of certain transition metal atoms^'^ 
or doubtful equilibrium geometries of Carbon clusters^^, 
and general, like Van der Waals forces not being pre- 
dicted. 

It has long been recognized that important insight can 
be gained by studying the asymptotic relationship be- 
tween correlation structure and one-body density in scal- 
ing limits^'^'^^. In this paper we focus on the strongly 
correlated limit of the exact Hohenberg-Kohn functional 
first investigated in Ref. 19-21 in which electron repulsion 
dominates over kinetic energy (yielding a natural coun- 
terpart to the Kohn-Sham kinetic energy functional^^). 

The resulting limit, which can be interpreted as an op- 
timal transport problem^'^, is still unwieldy from a com- 
putational point of view, since it requires the computa- 
tion of the full A^-point density of an A^-electron system, 
a function on R^^ . Here we give a simple but we believe 
fruitful reformulation as a minimization problem over 2- 
point densities subject to a representability constraint. 
This is similar to the well known formulation of the full 
quantum A^-body problem via representable 2-point den- 
sity matrices*'^^, but an important difference is that here 
it is only required that the 2-point density arises from 
a symmetric TV-point density rather than from an anti- 



symmetric, spin-dependent A^-particle wavefunction. We 
therefore speak of N-density representability. Density 
representability no longer mirrors the fermionic nature of 
electrons, refiecting the fact that a "semi-classical" limit 
has been taken of the Hohenberg-Kohn functional. 

We then establish a natural hierarchy of necessary rep- 
resentability conditions, and investigate the accuracy of 
the resulting reduced models as compared to the full 
strongly correlated limit. We focus on two test cases: 
first, a simple but illuminating 2-site, A^-particle model 
in which all representability conditions can be computed 
explicitly; and second, ab initio as well as self-consistent 
densities for the atoms He, Li, Be. A tentative conclusion 
is that the low order representability conditions already 
capture the main part of the correlations, at significantly 
reduced computational cost. 



II. STRONGLY CORRELATED LIMIT OF 
THE HOHENBERG-KOHN FUNCTIONAL 

The following counterpart to the Kohn-Sham kinetic 
energy functional was introduced in Ref. 19-21: 



Kt'^^W 



inf (*|T4c 



I* 



(1) 



Here the minimization is over electronic wavefunctions 
"if = ^(xi, si; . . . , xtv, Siv) which depend on N space 
and spin coordinates and belong to the usual space An 
of square- integrable antisymmetric normalized wavefunc- 
tions with square-integrable gradient, and the notation 
^ i-> p means that ^ has single-particle density p. The 
acronym SCE stands for strictly correlated electrons^°. 
Here we briefiy review known theoretical properties of 
1/,|^^ and previous approaches to compute it numeri- 
cally. 

Alternative constructions which plausibly yield the 
same functional are: 



i) semiclassical limit of the Hohenberg-Kohn functional: 



lim min ( *|r T + V^d* ), (2) 

ii) minimization over spinless bosonic wavefunctions $: 

(3) 
iii) minimization over A^-point probability measures: 



inf $|K=e|$/, 



mm 

PN^P, Pive'P^'" Jr3I 



Vr,, 



Pn- 



(4) 



Here Bn denotes the analogue of the space An for spin- 
less symmetric (bosonic) wavefunctions, P^™ stands for 
the set of symmetric probability measures on M^^, and 
Vee is the multiplication operator with the interaction 
potential 

yee(xi,...,XAr) = ^ Vco(Xi,Xj), (5) 

l<i<j<N 

where Wcc(x, y) = |x — y| . Formulae (2), (3) appear in 
Ref. 20, and (2), (4) are imphcit in Ref. 19. 

The minimum value in (2) with h = 1 is the exact 
Hohenberg-Kohn functional F^^ in the Levy-Lieb con- 
strained search formulation; so expression (2) is the semi- 
classical limit of the Hohenberg-Kohn functional. Note 
that the Kohn-Sham kinetic energy functional T^^ is 
obtained from F^^ by instead retaining the kinetic en- 
ergy operator T and neglecting the interaction term V^c 
Expression (3) is related to (1) by neglecting antisym- 
metry and spin, and to expression (4) by first noting 
that ($|Vj;o|$) = /Kie 1^1 and then replacing squares 
of spinless symmetric wavefunctions by their mathemat- 
ical "closure" , symmetric probability measures. 

Equality between the four expressions (1) - (4) was 
conjectured in Ref. 20 and has recently been justified 
mathematically^ . 

As noticed in Ref. 3 and 6, the last expression, (4), has 
the form of an optimal transport problem. In the standard 
setting of such problems^'^ originating from economics, 
one has N = 2, p2(x,y) corresponds to the amount of 
"mass" transported from x to y, T4e(x,y) is the "cost" 
of this transport, the one-body densities of x and y would 
be different from each other but prescribed a priori, i.e., 
/P2(x,y)dy = /9'^(x) and /p2(x,y)dx = p^{y), and 
minimization of J V^c p2 amounts to finding the most eco- 
nomical way of transporting the pile of mass p to p^ . 
In economics, the cost would typically increase rather 
than decrease with distance, prototypical examples be- 
ing |x-y| or |x-y|^. 

An interesting feature of minimizers is that they typi- 
cally concentrate on lower-dimensional sets (see Fig. 1). 
For N = 2, these sets have the form y = T(x). Phys- 
ically, this reflects the fact that given the position of 
the first electron, the position of the second electron be- 
comes deterministic in the strongly correlated limit (1). 



When p is radially symmetric, T is known explicitly in 
terms of the inverse of the radial distribution function 



i? K> 47r /q r'^p{r) dr. 



3,6,19-21 




FIG. 1. Optimal pair density p2 of Eq. (4) evaluated 
on (a;i, 0,0,yi, 0, 0) (green) for one-body ground-state den- 
sity p of helium (brown). The height of the green surface 
(1 + |VT(x)| )~^'^p(x) indicates the prefactor of the Haus- 
dorff measure on the set y = T(x). 

The minimization problem (4), and indeed any optimal 
transport problem, has two alternative formulations. To 
obtain the so-called Monge formulation, one makes the 
ansatz 



PAf(xi 



. ..,XAr) 
P(xi) 



N 



<5(x2-T2(xi))---,5(xjv-Tw(xi)) (6) 



for transport maps or co-motion functions T^ : M"^ —> M^ 
which preserves the one-body density p, that is to say 
It(a) P ~ IaP ^'^^ general subsets A C M"^ (see Ref. 6, 20, 
and 21 for physical and mathematical justifications). In 
fact, the ansatz (6) is not in general symmetric, so strictly 
speaking we should minimize over the symmetrizations of 
measures of form (6), but dropping the symmetrization 
does not alter the minimum value in (4). Or one passes 
to the so-called Kantorovich dual formulation 



Kr"[p] = 



sup 

, Yl «(xi)<Vco(xl,...,Xi^ 



pu (7) 



(see Ref. 3 for a mathematical justification and Ref. 16 
for a numerical scheme). 

The Monge formulation amounts to a spectacular di- 
mension reduction, in that the unknowns are A'^ maps on 
M.^ instead of one function on M.^'^ . Thus, when discretiz- 
ing M'^ by K gridpoints one has K ■ 3N instead of K'^ 
computational degrees of freedom. However, for A^ > 2 
it is not clear if the (symmetrized) Monge formalism cap- 
tures all minimizers of (4). See Section VIII for a coun- 
terexample when the Coulomb repulsion is replaced by a 
repulsive harmonic interaction. Moreover, previous nu- 
merical (and analytical) computation of V^o^^ using (6) 
is hitherto restricted to spherically symmetric densities 
or ID systems. This is because one has to deal with the 



infinite-dimensional nonlinear constraint that the T^ pre- 
serve p, and because the T^ are expected to jump along 
surfaces; in the radial case, the surfaces are believed to 
be concentric spheres^". 

The Kantorovich dual formulation, which has been suc- 
cessfully applied to non-spherical problems^^, cures the 
high storage complexity of the original formulation (4), 
but the inequality constraint in (7) that needs to be sat- 
isfied by u is still high-dimensional. 

It is then of interest to explore alternative ways of 
reducing the dimensionality of (1). The remainder of 
this paper is devoted to developing such an alternative 
approach, based on minimization over 2-point densities 
satisfying representability constraints. 



III. N-DENSITY REPRESENTABILITY 
AND REDUCED DENSITY MODELS 

We now derive a simple but we believe fruitful refor- 
mulation of the minimization problem (1). We begin by 
formalizing the notion of reduced densities. Throughout 
this section it is useful to work with the convention that 
all densities and reduced densities integrate to 1. We 
denote fc-point reduced densities with this normalization 
by Pk, to distinguish them from the customary fc-point 
reduced densities pk which integrate to the number of k- 
tuples in the system. Thus, given a symmetric iV-point 
probability density pn on M'^^, iV > 2, we define the 
associated one- and two-point reduced densities (known 
in probability theory under the name marginal densities) 
by 



Pi(xi) 

P2(X1,X2) 



R3(N-2) 



PAr(xi,. .. ,XAr) dx2- •• dxjv, (8) 
PAr(xi,. .. ,XAr) dx3- •• dxAT. (9) 



In particular, pi is related to the customary one-body 
density p by the formula pi — p/N . 

Typical iV-point densities occurring in the SCE limit 
concentrate on lower dimensional subsets (see Figure 1). 
Mathematically this does not pose real difficulties (and 
is a higher-dimensional analogue of the familiar charge 
distributions on surfaces in electrostatics). It just means 
that these densities should properly be regarded as prob- 
ability measures on K^^, not functions, i.e., they are not 
specified by pointwise values but by their integrals over 
sets. In the more general setting of probability measures, 
(8), (9) have to be replaced by 



dpi 



AxR3(«-i) 



dp 



N, 



(10) 
(11) 



dp2 = / dp AT, 

I A-X.B JAx_BxR3(w-2) 

for any subsets A,B C M.^. 

We employ the usual notation p^ '-^ Pi and ppf i-^ p2 
for the validity of Eq. (8) respectively (9). 



Definition III.l. (N -density representability) Let N > 
2. A probability density (or probability measure) P2 on M.^ 
is called iV-density-representable if there exists a sym- 
metric probability density (or probability measure) pj^ on 
M.'^^ such that p^ i—)- P2 • 

Examples 1) It is clear that p2 is 2-density representable 
if and only if it is symmetric. 

2) Any statistically independent measure p2(x, y) = 
Pi(x)pi(y) is TV-density representable for all N, since 
it is represented by the iV-body probability measure 
Pi(xi)---pi(x7v). 

3) The totally anticorrelated probability measure 

P2(x,y) = i(<5(x-A)J(y-S) 



-5(K-B)5{y-A)\, (12) 



A,Bg E'^, a y^ B is 2-density representable, but not 
3-density representable. That is to say, even though it is 
a symmetric probability measure on M^, there does not 
exist any symmetric probability measure p^ on M^ such 
that /p3(x,y,z)dz = p2(x, y). The reason is explained 
in Section IV. 

4) The previous example can be turned into a smooth 
one (see Figure 2). The smooth pair density 



P2(x,y) 



^(^(x 



A)ipiy-B) + ip{x-BMy-A)^ 

(13) 

with ip any nonnegativc function on M"' with J ip = 1 
and tp{z) — when \z\ > \A — B\ /2, is not 3-density 
representable, as we will show in Section V. 



P2^>^y) 




FIG. 2. Pair densities which are not 3-density-representable, 
such as the one depicted here (Eq. (13)), can be quite innocent 
looking. For further discussion of this example see Section V. 

The above definition immediately implies the following 
theorem. 



Theorem III. 2. Let N > M > 2. If a probability 
density (or probability measure) p2 on M^ is N -density- 
representable, then it is also M -density-representable. 

In other words, A^-density representability becomes a 
more and more stringent condition as N increases. 

Proof. If pj\j is a symmetric iV-body density which rep- 
resents p2 , then 



p(xi,...,xm) 

= / PAr(xi,... ,XA./,...,XAr)dX7\/+i • •• dx 

JmHN-M) 



N 



is a symmetric M-body density which also represents p2- 

D 

With the help of the concept of density representabil- 
ity, we can exploit the fact that the Coulomb potential 
Vee in (4) only involves pair interactions to reformulate 
the many-body optimal transport definition (4) of V^J~^^ 
as a standard (two-body) optimal transport problem with 
a constraint. This result does not depend on the Coulom- 
bic form of the interaction potential Wcc- 

Theorem III. 3. (SCE energy via density representabil- 
ity) For any given single-particle density p of an N- 
electron system, 



Kt""[p] 



mm 

P2 N -density-rep. 



N 



VccP2, 



(14) 



where Woo is any pair potential which is symmetric (i.e., 
w(x, y) — v(y,x)), and the minimization is over proba- 
bility densities P2 on M^ . 

Here the appearance of the normalization constants N 



and 



/N\ 



is due to the fact that p integrates to N, not 1, 



whereas p2 integrates to 1, not 



(N\ 



Proof The proof is similar to the famous proof of Levy 
of the Hohenberg-Kohn theorem. For any symmetric Ap- 
point density pN with p^ ^-^ P2, we clearly have 



Vcc Pn 



VcoP2\ 



(15) 



that is to say the electron-electron energy only depends 
on the two-body reduced density of pjv- We can therefore 
usefully partition the minimization in (4) into a double 
minimization, first over pjv subject to fixed P2, then over 
P2- 



mm 

Pm^p/N 



VccPN 



mm mm 

P2i-^p/N, Pn>->P2 

P2 A^-dcnsity-rcp. 



mm 

P2i-^p/N, 

P2 A^-dcnsity-rcp. 



VccPN 



VccP2, 



the last equality being due to Eq. (15). 



D 



The formula for Vcc*^^ in Theorem III. 3 together with 
the necessary conditions for iV-density representability in 
Theorem III. 2 immediately suggests a natural hierarchy 
of approximations. For any given single-particle density 
p of an TV-electron system, let us define 



K 



SCE,*: 



[P] 



mm 

P2<-^p/N, 
P2 fc-dcnsity-rcp. 



N 



VccP2, 



fc = 2,3. 



(16) 



That is, we replace the requirement that p2 be A^- 
representable by the weaker requirement that it be k- 
representable for some k < N. This enlarges the set of 
admissible p2^s in the minimization, leading to the fol- 
lowing chain of inequalities 



ySCE,: 



< 



V' 



SCE.fc 



[P] 



< 



min / Wcc P2 

P2^^p/N 



VcT^^^ip] 

II 



(17) 

We call V^J~^^'''' the order-fc approximation of the SCE 
energy. The lowest-order approximation V^^ ^'^ corre- 
sponds to solving a classical (two-body) optimal trans- 
port problem with Coulomb cost (yielding the functional 
introduced in Ref. 6 which we called F'~^'^[p]), whereas the 
order-N approximation 14|^^'^ recovers the exact SCE 
energy. Physically, the intermediate functionals V^c^^''^ 
can be thought of as reduced models for the energy of 
strongly correlated electrons which take into account k- 
body correlations. 

The unknown in the order-A; approximation is a /c-body 
density, so the computational cost increases steeply with 
k; e.g., discretizing each copy of M.'^ by K gridpoints 
leads to a iiT'^-point discretization for M^'^. The practical 
value of our reduced models therefore depends strongly 
on whether low-order approximations are already capa- 
ble of capturing the main part of the full SCE energy. A 
tentative answer is that they are, as we will document in 
the next two sections. 

Theoretically the parameter k in (16) can also be chosen 
bigger than N, a particularly interesting question being 
what happens when fc — >■ oo. In the present paper, we 
will only answer this question for the model densities (18) 
below. A general discussion will appear elsewhere. 

Finally, we remark that density representability of a pair 
density is obviously a necessary condition for the familiar 
(wavefunction) representability of any two-body density 
matrix which gives rise to this pair density. Analyzing the 
relationship between this necessary condition and com- 
mon representability conditions from density matrix the- 
ory such as the P, Q and G conditions'*'^''' lies beyond 
the scope of this paper. 



IV. MODEL PROBLEM: N PARTICLES 
OCCUPYING 2 SITES 

In this section we analyze a model system in which 
the particle positions are restricted to 2 sites, to gain ba- 
sic insights into what it means for a pair density to be 
/c-density representable and into how the resulting func- 
tional (16) depend on k. The single-particle density of 
such a system has the form 



•i8)^i«, «i 



, zat G {A, B}. The set of symmetric 



N 



{l~t)6A + tSB, 0<t<l, 



(18) 



where N is the number of particles and A and B are 
two different points in M.^. This model density, while of 
course very simplistic, can be regarded as a toy model for 
the electron density of a diatomic system in the regime 
when the interatomic distance is much larger than the 
atomic radii. If we don't want to allow fractional occupa- 
tion numbers of the sites, t would be restricted to integer 
multiples of 1/N, but since this makes little difference 
to the analysis, we might as well allow real occupation 
numbers. 

Our first goal is to compute explicitly the set of A^- 
representable 2-point probability measures for our 2-site 
system. The TV-point probability measures on M.^^ whose 
single-particle density has the form (18) for some t are 
the measures of form 



7=(zi,...,i„)e{A,S}" 



(19) 



with a/ > 0, J2i '^i — 1' ^H-d correspond to the 
probability measures on the discrete 2-site, A^-particle 
state space {A,B}'^. We use the following notation 
for the different sets of probability measures of interest: 
V{{A,B}'^) denotes the set of probability measures on 
{A,B}^, i.e., all measures of form (19); r''y"'{{A,B}'^) 
is the set of such measures which are symmetric, i.e., 
'^(ji,...,ijv) is ^ symmetric function of its arguments 
(ii, . . . , iM); and r^-''''Pi{A, B}2) stands for the set of A^- 
density-representable probability measures on the two- 
body state space {A,B}^, i.e., those probability mea- 
sures on {A, B} which arise as marginals (9) of some 
PN e 7'"y"({A,B}^). In particular, the 2-density- 
representable probability measures are those measures of 
form 

P2 = aAASA 1^ 6a + ubbSb ® Sb 

+ aAB^A <E)Sb + OiBA^B ® 5a (20) 

which satisfy the trivial conditions of nonnegativity, nor- 
malization, and symmetry, 

aij > for all i, j, ^ aij = 1, uab = asA- (21) 
hi 

It is clear from the explicit representation (19) that 
'P{{A,B}'^) is the convex hull of its extreme points 



iV-point probability densities satisfies 'P'''""({A, B}^) = 
Sis!'P{{A,B}^), where Sn is the symmetrizer 

{SnPn){AiX----kAn) = — ^pjv(^a(l)X---xA^(jv)) 



and the sum runs over all permutations. It follows that 
T^^y^dA, B}^) is the convex hull of the elements SN5i^ ® 
• • • ® (5,„, and that P^-"p({A, BY) is the convex hull of 
their two-point densities, 

'P^"'°P({A, B}) = convex hull of the measures 



\P2 : «!,■ 



, ITV 



e{AB}}, (22) 



where here and below, P2 " denotes the two-particle den- 
sity oi Pn- To compute these two-point densities, we use 
an averaging formula which can be shown by an elemen- 
tary computation: first symmetrizing and then taking 
the two-point density is the same as taking the average 
over all two-point densities, 



(SnPn) 
P2 



1 

\2) l<i<j<N 



E 



Pn dx, 



I], 



(23) 



where x^- denotes the list of coordinates Xi , . . . , xjv with 
Xi and Xj omitted, and pn is any A^-point probability 
measure. Now consider a measure of form pN = Si^ (g) 
■ ■ ■ (S Si^, and let K = ]),{ij \ ij = B}, i.e. the occupation 
number of site B. Note < K < N. By the averaging 
formula (23), and using the abbreviated notation 6^ — 



P2 



— P2 



'<»S^f 



Q 



r /TV - K 



5a® 5a 



K 



K{N -K) 



{5a ® 5i 



(24) 



Note that the resulting two-point marginal does not de- 
pend on the ij, but only on the occupation number 
K e {0, ...,Af}. Equations (22), (24) give the follow- 
ing final result: 

Theorem IV. 1. The set of N -representable 2-point 
measures, V^~^''^{{A, B}), is the convex hull of the K-\-l 
measures given by the right hand side of (24), where K 
runs from to N . 

This set is plotted in Figure 3, for different values of 
N. 

Next we show that, as suggested by Figure 3, when N 
gets large the iV -I- 1 extremal pair densities in (24) ap- 
proach a certain very interesting continuous curve. To 
see this, let us re-write formula (24) in terms of the nor- 
malized occupation number t = K/N S [0, 1] instead of 



, + Q 




J^AA 



FIG. 3. The set of A''-representable pair densities of form 

aAASA^SA+CtABSA^SB+CtBASB^SA+aBBSB^SB ioi N = 2 

(red), TV = 3 (green), N = 6 (blue), and iV = 20 (light blue). 
The coefficient oab + ctBA on the vertical axis indicates the 
weight of the anticorrelated contribution (12); its maximum 
representable value decreases with A'^. When N = 2, only the 
trivial conditions (21) are present. Remarkably, as A'^ gets 
large the upper boundary of the representable set approaches 
the curve given by the mean field densities, i.e., P2 ~ Pi ^ Pi 
for some pi (grey curve), see Eq. (25), (27). 



K, and separate the coefRcients into iV-independent and 
lower order terms. An elementary calculation shows that 



K - 1 



t 



1-t [N - K)-l 



N -I 7V-1' N -1 

and, abbreviating 6i <Si Sj by dij , 



(1-t) 



t 



N -I' 



P2 



(1 - tfSAA + t^SsB + t(l - t){SAB + 5ba) 



-Pt' 



N -1 



i-s^ 



AA - ObB -I- OAB -I- Ob A 



(25) 



But the first term is precisely the mean field approxima- 
tion to the pair density of the state SnS^^^5§ obtained 
from its single-particle density 

P^f"^ '''^^{l-t)SA+t6B, (26) 

namely 

P^^ = ((1 - t)SA + tSs) ® ((1 - t)SA + tSs) . (27) 

The second term in (25) is a correlation correction which 
depletes the "ionic" terms 5a ® Sa and 5b ® 5b in favour 



of the "anticorrelated" terms 5a ® 5b and 5b® 5a- This 
correction is large for small TV (and even completely re- 
moves the ionic terms when N — 2 and t — 1/2), but 
vanishes in the limit TV — > cx) at fixed occupation number 
t. 

In particular, we have established the following 

Theorem IV. 2. A pair density of form (20) is N- 
representable for all TV if and only if it lies in the convex 
hull of the mean field densities, or ~ by inspection of Fig- 
ure 3 - if and only if it is a convex combination of two 
mean field densities. 

The "primal" description of TV-representable pair den- 
sities as the convex hull of explicit extreme points can be 
easily turned into an equivalent "dual" description via 
inequalities. We only give the result in the cases A^ == 3 

and TV = oo. 

Corollary IV. 3. A pair density of form (20) is 3- 
representable if and only if it satisfies (21) and the linear 
inequality 



UAB + OIBA < 2(aAA + "Bs), 



(28) 



and N -representable for all TV if and only if it satisfies 
(21) and the nonlinear inequality 



aAB + a.BA < 2y/aAA ■ aBB- 



(29) 



To derive (29), one first shows that uab + oiBA < 
2{aAA + oiAB){oiBB + olba)- Thanks to Eq. (21) this 
is a quadratic inequality for olab + olba and solving it 
yields (29). 

The physical meaning of Eq. (28) is that at most 2/3 
of the mass of p2 can sit on the non-ionic configurations 
{A,B) and {B,A). The meaning of Eq. (29) is that the 
total size of the non-ionic contributions cannot exceed its 
size in the mean field pair density formed from its single- 
particle density. 

The above results can be used to determine the hierar- 
chy of approximate TV-particle functionals V^,^^^'*^ intro- 
duced in (16) on densities of the form (18). In fact, the 
exact Coulomb potential Wcc(x, y) = |x — y| no longer 
makes sense in the context of these model densities since 
multiply occupied sites would lead to infinite energy, so 
we replace it by an appropriately regularized interaction, 
with the property that 



M.A) 



> Uab 



^diag 
= Woc(A,S) 



v,,{B,A). (30) 



Here t/diag and Uab are effective parameters for same-site 
and different-site repulsion, and the inequality C/diag > 
Uab preserves the repulsive effect that the interaction 
potential decreases with interparticle distance. Hence 
the two-point densities p2 which compete in the varia- 
tional definition (16) prefer the different-site configura- 
tions (x, y) = {A, B) and (x, y) = (i?. A) over the ionic 



configurations {A, A) and {B,B). Consequently the op- 
timizing p2's with one-point density (18) are those k- 
representable 2-point densities of form (20) which have 
one-point density pt (this fixes their position in direction 
of the baseline in Figure 3, because t = \{aBB — ctAA)) 
which maximize the coefficient uab + c^ba, i-e., lie on 
the upper boundary of the representable set in Figure 3. 
When t is an integer multiple of 1/fc, i.e., t = K/k, 
K ~ 0,1,. . . ,k, the optimizing p2 is thus precisely given 
by formula (25) with N replaced by k. It follows that, 
denoting the right hand side of (18) by pt. 






(^] (c/diag • [t^ + (1 - tf] + Uab ■ 2t{l - t)) 

2 j (C/diag -C/ab)^^^, 

t^K/k, K = 0,l,...,k. (31) 



1 T/SCE,kr^ 1 
(^^ee [Pt\ 




Pt 



For intermediate occupation numbers t with t_ = {K 



l)/k < t < K/k = t 



+ ^ 



K 



1, 2, . . . , fc, the upper 



boundary of the fc-representable set is given by the lin- 
ear interpolation between the p2's coming from t_ and 
t+, and hence so is the resulting value of V^^^'^ ■ The 
weights of the contributions from t± are the same as 
the interpolation weights for the single-particle density, 
Pt^{K~ kt)pt_ + {kt - {K - l))/9t+, and so 



FIG. 4. Density-representability approximation of order k to 
K:o , on densities of form pt/N = {1 — t)5A -\- tSs (various 
values of k). The Coulomb interaction has been replaced by 
the regularized interaction (30). Red, green, blue, and pink 
corresponds to fc = 2,4,6, 11. The piecewise linear structure 
is an exact feature of the results (31), (32). The order-fc 
approximation equals the exact Kc for k = N particles. 
Light blue curve: mean field energy (see text). 



K 



SCE,fc 



N = 



(if-A;t)Kf^''=b( 



(kt - [K 



K-l)/k 

i))Kf^''= 



K 



- 1 K 



[PK/k], 

K = 1. 



,fc. (32) 



The reduced SCE energies (31), (32) are plotted in Fig- 
ure 4, for different values of k. 

Finally let us calculate and physically interpret the 
large-fc limit. It is clear from our explicit results that the 
limit is just given by the first part of the right hand side 
of Eq. (31). This part is nothing but a (Hartree-type) 
mean field energy. 



iimv;fE.'=N = r 

/v — rCXD \ ^ 



Pt ^ Pt 



= 1 



J[P] 



N N 
«cc(x,y)p(x)p(y)dxdy 



J[Pt 



(33) 



Here the prefactor 1 — 1/iV is a self-interaction correc- 
tion, i.e., the approximation via density representability 
of infinite order remembers that there are only ( 2 ) in- 
teraction terms, not N'^ /2. In other words, remarkably, 
the infinite-order approximation to the SCE functional is 
nothing but the self-interaction-corrected mean field en- 
ergy, even though mean field approximations played no 
role in the construction of the reduced SCE functionals 
(16). 



After completing the above (elementary) analysis, we 
learned that results similar to, and in fact much more 
general than, Theorems IV. 1 and IV. 2 are well known in 
probability theory, more precisely in the theory of "ex- 
changeable sequences" of random variables^ '^. This the- 
ory, which appears to be hitherto disconnected from DFT 
(as well as wavefunction representability), entails a classi- 
fication going back to de Finetti^ of symmetric probabil- 
ity densities Poo(xi, X2, . . . ) in infinitely many variables. 



V. NECESSARY CONDITIONS FOR 
DENSITY-REPRESENTABILITY 

The results on two-state systems in the previous sec- 
tion immediately yield necessary conditions on iV-density 
representability for general pair densities p2 on M^. To 
this end, let us introduce the following integrals of p2 
associated with any partitioning of M.^ into two disjoint 
subsets Qa E^nd fi^: 



P2, i,je{A,B}. 



(34) 



fli X Uj 



Theorem V.l. Let p2 be any pair density (or measure) 
on M.^ , normalized so that J p2 = 1- If P2 is N-density- 
representable, and Ha, ^b is any partitioning ofM.^ into 
two subsets, then the associated 2-site pair density 



OtAA^AA + (XAB^AB + CUbaSbA + O^BbSbB, 



(35) 



with aij as in (34), is also N -representable, that is to say 
it belongs to the set 'P^~^'^'''{{A, B}) computed explicitly 
in Theorem IV. 1 and depicted in Figure 3. 

Proof. If pn is an A^-point density (or measure) on M?-^ 
which represents p2, then the associated 2-site, A'^-point 
density 



ii,...,iN&{A,B} 



•iiv "ii 



with a. 
represents the 2-site pair density (35) 



11---IN - I Pn, (36) 



n 



Example The smooth anticorrelated pair density (13) 
in Example 4 of Section III (see Figure 2) is not 3- 
representable: choose ^Ia, ^b to be the half-spaces of M'^ 
whose boundary bisects the line segment from A to B, 
that is to say rj^ = {x € R^ | x • (B - A) < i\f • (S - ^)}, 
Ob = {x e E3 I X • (S - A) > Af • (B - A)}, where 
M ^ {A + B)/2. By construction, in this case aAA = 
<xbb — 0, aAB ~ OLBA = 1, so the associated 2-site pair 
density (35) is not 3-density representable, as shown in 
the previous section (see Figure 3). 

In fact, for iV = 3, the dual description (28) of the A^- 
representable pair densities for the 2-site system yields 
the necessary condition of Theorem V.l directly in the 
form of the following inequality: 



P2^ 



flAXfl^ 



P2<2 



ObX^a 



P2 



QAXfi^ 



P2 



Ub 'xilf 



(37) 

Physically, this condition says that the total probability 
to find a particle pair in the "anticorrelated" regions VI a x 
r^B and VLb x VIa can be at most twice as large as the 
probability to find a pair in the "ionic" regions Q.a x ^a 
and ^B X ^B- 



VI. BEHAVIOUR OF THE REDUCED 
MODELS ON AB INITIO DENSITIES FOR 
SMALL ATOMS 

We now investigate the effect of the hierarchy of rep- 
resentability conditions (Section III) for the atoms He, 
Li and Be. The ab initio single electron density p is ob- 
tained from a full configuration interaction (FCI) calcu- 
lation with Slater- type orbitals (STOs)^". The approxi- 
mate interaction energy V^^^''^ of an TV-electron system 
can be obtained by simulating a fictitious fc-electron sys- 
tem: directly from Eq. (16), we have 



K 



SCE./c 



[P] = 



(I) 
(2) 



V' 



SCE 



^P] 



(38) 



The energy on the right hand side is obtained by the same 
method as in Ref. 20; in particular, the jump surfaces of 
the maps T^ in (6) are assumed to be concentric spheres. 





He Li 


Be 


fc = 2 


0.711906 1.17881 


2.21361 


fc = 3 


1.48426 


2.88229 


fc = 4 




3.24853 


exact 


0.954988 2.30755 


4.51366 


J 


2.67842 4.11866 


7.30589 



TABLE I. The calculated values of V^^'^^''' using the fc-density 
representability approximation (see also Fig. 5). Each diago- 
nal entry is the true K:c • An "exact" (FCI) value and the 
Hartree term J are also shown, for comparison. 



Fig. 5 compares the fc-density representability approx- 
imation V^^^'^[p\ with Ko^^[p], with the "exact" value 
from the FCI calculation, and the Hartree term J. V^^*^^ 
underestimates the exact value, whereas J overestimates 
it. The order-fc approximation is already reasonably close 
to Vp^*^^. The corresponding numerical values are sum- 
marized in Table I. 
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FIG. 5. Vel'^^'^H obtained by the fc-density representabil- 
ity approximation (blue) and the true Voc [p] (cyan). The 
green curve shows the "exact" {>I'|Kc|'I') within an STO FCI 
ansatz space. The mean-field Hartree term J (orange dashed) 
overestimates the exact value, as expected. See Tab. I for the 
numerical values. 



VII. A SELF-CONSISTENT KOHN-SHAM 
COMPUTATION COMPARING EXACT AND 
REDUCED SCE 

The SCE formalism has the potential to become 
an important ingredient in the design of exchange- 
correlation functionals for strongly correlated electron 
systems. Thus we investigate the SCE approach in the 
context of a Kohn-Sham self-consistent field calculation 
for atoms with the total energy functional 



E[p]=Tks[p\ + V' 



SCE.fe 



[P] 



t(x)p(x)dx. (39) 



Here Tks is the Kohn-Sham kinetic energy functional 
and Uoxt is the external nuclear potential. Previous self- 
consistent field calculations with the SCE functional were 



carried out in Ref. 13 and 16 for a ID quantum wire, 
where in the weak confinement regime the SCE func- 
tional becomes asymptotically exact. For 3D atomic sys- 
tems considered here, replacing J -I- E^c by Kc'"^''' (or 
even by the exact SCE functional T4o^^) presumably does 
not yield physically accurate results due to the missing 
influence of kinetic energy on p2, but our calculations 
illustrate the effect of the fc-density approximation. 

The Kohn-Sham equations define a nonlinear eigen- 
value problem 



H[p]ipi = s.iip. 



(40) 



p(x)=^|^,(x)^ /'^,(x)>,(x)dx = <5,„ (41) 

i=l ■' 

where the Hamiltonian H[p] itself depends on the den- 
sity p. For an atom with nuclear charge Z, and the den- 
sity functional (39) with k — N (exact SCE), the single- 
particle Hamiltonian (in atomic units) reads 



H[p] = -lA-^+u[py 



(42) 



The term — |A — Z/\x.\ is the hydrogcn-likc single- 
particle Hamiltonian, and u[p] is the Kantorovich poten- 
tial, i.e., the maximizer of (7), which enters because for- 
mally, {SV^J~^^/6p)[p] = u[p\. Note that changing the po- 
tential in (42) by an additive constant would not change 
the Kohn-Sham orbitals, but choosing precisely u[p] has 
the virtue that the Kohn-Sham eigenvalues sum to the 
system energy E[p], as is easily inferred from (40), (7). 

The Kantorovich potential agrees up to an addi- 
tive constant with the effective SCE potential v[p] con- 
structed in Ref. 20. The latter can be defined by^" 

Vz;[p](x) ^-Y: ," ^fi . , lim v[p]{^) ^ 0. 

,^2|X-T,(x)| IxKoo 

(43) 
The Ti are the transport maps in Eq. (6), which deter- 
mine the positions of the remaining electrons given the 
position of the first electron, and solely depend on the 
density p. The additive constant is now easily obtained 
from (7): u[p] = v[p] + C with 



C 



E 



1 



^ fe|T.W-T,(x) 



■dx— / /9(x)w[/9](x) dx. 
(44) 



Following Ref. 20 we assume on physical grounds that 

(45) 



iV- 1 
w[p](x) ^ — I — I — as |x| — > cxo. 



even though we do not know of a mathematical proof. For 
charge-neutral atoms with N ^ Z, the Hamiltonian (42) 
can be re-written as 



^[^]— ^R+^- 



u[p](x)- 



N-1 





He Li Be 


fc = 2 


-2.74058 -5.24439 -8.12061 


fc = 3 


-5.10739 -7.90984 


fc = 4 


-7.79889 



(46) 



TABLE H. Kohn-Sham energy (sum of eigenvalues) obtained 
by a self-consistent field iteration with the SCE Kantorovich 
potential and Hamiltonian in Eq. (46). 



such that the last term is expected to decay faster than 
1/ |x| due to the asymptotic relation (45). 

Fig. 6 shows the self-consistent densities and corre- 
sponding SCE potential v of helium, lithium and beryl- 
lium, for the exact SCE potential as well as its fc-density 
approximation (obtained from simulating a fictitious k- 
electron system, see Section VI) . All densities are normal- 
ized to N. For each fc, v[p] is rescaled by ^^3f to match 
the asymptotic expansion (45). Table II summarizes the 
numerical Kohn-Sham energy (sum of Kohn-Sham eigen- 
values, with doubly occupied orbitals due to spin). Note 
that the fc-density representability approximation of the 
energy is below the true value and increases with fc. 



VIII. EXAMPLE OF A MINIMIZING 
TV-POINT DENSITY NOT OF SCE FORM 

Since much of the theoretical and numerical work on 
the minimization problem (4) relies on the (plausible but 
nontrivial) ansatz (6), it is of interest to understand its 
precise status with respect to minimization over arbitrary 
A'^-point probability measures. 

For N — 2 the ansatz is known to be exact, in the 
sense that the minimizing 2-point probability measure 
is unique and of the form (6) (Ref. 6, following ear- 
lier work^^ in the optimal transportation literature on 
pair interactions Uco which increase with interparticle dis- 
tance) . 

For N > 2, physical arguments^*^ suggest that there 
should always exist a minimizing A^-point probability 
measure of this form, so in particular restricting the min- 
imization in (4) to Pat's of form (6) should always give 
the correct value of the functional V^l^^ [p] , but there is 
no rigorous proof of this conjecture. 

Finally, there is the question whether for A^ > 2 the 
ansatz (6) yields all solutions. Here we are not aware of 
any convincing arguments (be they physical or mathe- 
matical), one way or the other. The following counterex- 
ample demonstrates that for N > 2, (6) does not yield 
all solutions if the Coulomb interaction is replaced by a 
negative harmonic oscillator interaction. This is a new 
and somewhat surprising effect which only appears when 
N > 2; for N = 2 the negative harmonic interaction was 
already considered'^ '^° in connection with the SCE func- 
tional and the two interactions were shown to behave in 
exactly the same way. 

Our counterexample does not imply that the ansatz 
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FIG. 6. (a, c, e) Radial part of the self-consistent density (blue) for tlie lielium, lithium and beryllium atom with the 
SCE exchange-correlation functional _Exc = KI — J ■ The lighter blue curves correspond to the fc-density representabihty 
approximation, and are visually indiscernible from the (exact) k = N case, (b, d, f) The SCE (alias shifted Kantorovich) 



potential v[p] (blue) corresponding to the self-consistent density on the left, rescaled by 
asymptotic expansion — ^j^ in Eq. (45) . 



The gray dashed line shows the 



(6) does not capture all miiiim.izers of the true Coulombic 
SCE problem (4), but it means that if it does, this must 
be because of some special Coulombic features. 

The counterexample is best discussed in the context of 
recent work in the optimal transport literature on A^- 
body optimal transport problems in M.'^^ with a gen- 
eral nonnegative interaction potential or "cost function" 

Voo(xi,.. .,XAr), 



mm / VccPn, 

PjvU-P JjJdJV 



sures on 



^dAf ^g shown in Ref. 17, minimizers have to 
concentrate on subsets whose dimension is bounded in 
terms of the signatures (the number of positive, negative 
and zero eigenvalues) of certain symmetric matrices de- 
rived from the mixed second order partial derivatives of 
V^o- Let G be the off diagonal part of the Hessian of Voc. 
More explicitly, if 



D'Vo. 






where the minimization is over iV-point probability mea- denotes the d x d matrix of mixed second order partials 
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with respect to x^ € 



and Xj G 



G = 



^^,x.Kc 



D^ V 



-'^XjnXi '''^cc ^x„^x2 ^oe 



we have 



^x,x,„K=c 

D'^ V 







(47) 



Note that G is a block matrix; each entry in the preceding 
formula denotes a d x d block. Now, as a symmetric 
NdxNd matrix, G has Nd real eigenvalues, counted with 
multiplicities. Let A+,A_ and Aq denote, respectively, 
the number of positive, negative and zero eigenvalues of 
G at some point x ~ (xi,X2, . . .xjv) € M^'*; note that 
A_|_ + A_ + Aq = Nd. Then, near x, Theorem 2.3 in 
Ref. 17 implies that the support of minimizers (the subset 
on which they are nonzero) is contained in a subset of 
dimension Aq + A_. 

For the Coulomb interaction '^i^j\^i^^j\ , a 
straightforward calculation implies that 



1 



X,; X^ 



, yX.^ Xj ) (^X^ Xj ) 



where / is the d x d identity matrix. The signature of 
G, however, may change depending on the point x. One 
can show that (except at special points) d < Aq + A_ < 
(N — l)d, meaning that the dimension of the support of 
the solution can be no more than {N — l)d. In particular, 
for N = 2, this yields an alternative justification of the 
ansatz (6). 

While the preceding result is only an upper bound on 
the dimension, it is nevertheless a useful guideline for 
constructing high-dimensional minimizers, since the G- 
matrix of the cost must then necessarily have a large 
number of nonpositive eigenvalues. 

For ease of analysis of the G-matrix, consider now a 
cost of pair potential form, Eq. (5), with Vcc symmetric 
and quadratic. The dxd block D^^ Vee is then indepen- 
dent of i, j, and x, and the signatures of G can be com- 
puted explicitly. The maximum number of nonpositive 
eigenvalues occurs for the negative harmonic oscillator 
interaction 



'^cc l^Xj, Xj j 

In this case each D?:.^.Vcc 



(48) 



G = 



2/ 
2/ 

2/ 2/ 



21, and 

2/ 
2/ 







Now for any v E M'^, {v,v, . . . , v) is an eigenvector with 
eigenvalue 2{N — 1), while the vectors 

{v, -V, 0, ... , 0), (v, 0, -t;, 0, . . . , 0), . . . , (i;, 0, . . . , 0, -v) 



are eigenvectors with eigenvalue —2. This implies that 
A_ = {N — l)d, A+ = d and Aq = 0. Therefore, mini- 
mizers have at most {N — l)d-dimensional support. We 
now show that this bound is sharp for this cost function; 
that is, there actually are minimizers which are strictly 
positive on {N — l)(i-dimensional sets. 

Example Replace the Coulomb interaction with the 
negative harmonic oscillator interaction (5), (48). Let 
Pn be any symmetric measure on M.^^ which is con- 
centrated on the 3(A^ — l)-dimensional surface {xi -|- 
X2 -l- • • • -l- x^v — 0}. Then this measure is optimal 
for the corresponding single particle density p(xi) = 

^^/l[S.3(iV-l) /5w(xi, X2, . . . , Xjy) dx2 dXg, . . . dXAr. 

To see why, note that by a simple computation 



Vcc = - ^ |xi - : 



i<j 



1 



Xl 



Xat 






Hence for any p^ with one-body density p 

VccPn = / |xi H hxjv|^/7jv + - / |xi|^p(xi) dxi. 

The first term is minimized if and only if p^ is zero out- 
side the surface xi + • • • + x^v = 0, and the second term 
only depends on the one-body density p. Since p^ van- 
ishes outside this surface, it is a minimizer. 

The above example is in fact a special case of a result 
in Ref. 17. The interested reader is encouraged to consult 
17 for further results on the dimension of the support of 
optimizers. 

IX. CONCLUSIONS AND OUTLOOK 

We have reformulated the strongly correlated limit of 
density functional theory via 'W-density representabil- 
ity", i.e., the requirement that the pair density comes 
from a symmetric A^-point probability measure. This 
formulation gives rise to a natural hierarchy of approx- 
imate models, in which one relaxes this requirement to 
the existence of a representing symmetric fc-point density 
with k < N. In this paper we have presented a com- 
putational method for the approximate models which is 
akin to a wavefunction method, in that the representing 
fc-point density is resolved. One of the numerical find- 
ings we did not anticipate is the extreme robustness of 
self-consistent Kohn-Sham densities with respect to the 
fc-density approximation. 

For low fc, a promising route towards extending our 
methods to spherically asymmetric systems is the direct 
computation of the Kantorovich dual potential^^. In the 
future, if a more direct understanding of the main con- 
straints on p2 implied by fc-representability can be ob- 
tained, one could also envision a dual approach akin to 
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reduced density matrix methods in which one would solve 
a constrained linear programming problem for the pair 
density. 

Finally, another interesting issue raised by this work 
is to clarify the somewhat surprising connection between 
the SCE formalism and the mean field approximation 
suggested by our study of the two-site system in Sec- 
tion IV. 
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